Population Dynamics in the Penna Model 
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We build upon the recent steady-state Penna model solution, Phys. Rev. Lett. 89, 288103 (2002), 
to study the population dynamics within the Penna model. We show, that any perturbation to 
the population can be broken into a collection of modes each of which decay exponentially with its 
respective time constant. The long time behaviour of population is therefore likely to be dominated 
by the modes with the largest time constants. We confirm our analytical approach with simulation 
data. 

PACS numbers: 87.23.-n, 87.10. +c 



I. INTRODUCTION 

Genes, mutation, evolution and ageing have been top- 
ics of intensive research 0, H 0, El > particularly after 
the recent Genome project. In 1995, a bit-string com- 
puter simulation model for population evolution was in- 
troduced by Penna [|| which successfully encompassed 
all those elements. The Penna model, essentially a mu- 
tation accumulation model, has been so successful that it 
has rapidly established itself as a major model for popu- 
lation simulations 0] . Recently, analytical solutions have 
been presented for the steady states of the Penna model 
@ , shedding insights into the inter- relationships be- 
tween various Penna parameters of the simulation model. 
Here, we build upon this analytical framework to study 
populations undergoing steady growth or decline and to 
study the transient behaviour of the model when the sys- 
tem is away from its steady states. Previous attempts 
to analyse the population dynamics did not have the ad- 
vantage of the full analytic steady state solution @ , 
and consequently could not fully explore the dynamics of 
the model. We find that the fluctuations in a population 
away from steady state can be decomposed into a collec- 
tion of modes each of which decays exponentially with 
its respective time constant. The long time behaviour of 
a population will therefore be dominated by the modes 
with the largest time constants. Our analytical results 
are confirmed by comparison with simulation results. 



II. THE SIMPLE PENNA MODEL 

In the Penna model 6] , an organism's genome is repre- 
sented by a bit-string. Organisms age in time steps and 
at each timestep an organism reads the corresponding 
bit from the bit string (e.g. 2nd site at age 2). If a site 
contains a 1 the organism develops a disease, and once it 
has accumulated T diseases it dies. In any timestep an 
organism can reproduce with probability b. The child's 
bit-string is a copy of the parent's with a probability m 
of each site mutating into a 1. Positive mutations are 
rare in nature so a 1 mutating into a is forbidden in 
the model. Variants of the Penna model exist in which 
an organism can only reproduce between certain ages, 



and in which there is an external death rate giving each 
organism a genome-independent chance of death in any 
timestep. 

Traditionally ^(| the Penna model is implemented 
computationally and the population controlled by the 
use of a Verhulst factor [Tl| which controls either the 
birth rate or external death rate. The bit string is 32/64 
bits long to enable bitwise manipulation on integer type 
variables. The finite length bit string is an artifice of 
simulation and is not an important consideration when 
approaching the model analytically. A solution to the 
Penna model in the steady-state has been developed and 
is capable of dealing with a wide range of modifications 
to the model, namely arbitrary birth and survivability 
functions @- 

By building upon the steady-state solution it is pos- 
sible to consider dynamic behaviour. The simplest form 
the steady-state solution takes is for the simple Penna 
model in which there is no non-genetic source of death, 
an organism dies after a single disease (T = 1) and can 
reproduce with equal probability at any point during its 
life. For simplicity we present our dynamics analysis 
within this simple Penna model, since it is straightfor- 
ward to generalise to the case of T > 1 0. The steady- 
state solution to the simple Penna model is given in 
brief below: 

An organism within the Penna model can be uniquely 
characterised by its age x and the number of Os on its bit- 
string before the first 1. The number of Os determines 
how long the organism lives and is termed its string- 
length I. Where rij(x, I) is the number of organisms with 
age x and string-length I at time-step j 

oo 

n j+1 (0,l) = be- 0l ^rnj(x,l) (1) 

x=0 

oo oo 

l'>l x=0 

where = 1 — m is the probability of avoiding a mu- 
tation. In the stationary-state n J+ i(x, I) — rij(x,l) and 
given that an organism with string-length I lives for I 
timesteps, the sum over ages of n(x, I) is / x n(0, 1). The 
sum over ages is written as n(l) = n(x, I). Then eqn. 
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FIG. 1: Lifespan distribution for a simple Penna model with 
'max = 30, P — ^j. Analytical results (x) are compared with 
those from simulation (□ ). Simulation size 10 7 . 



||TJ simplifies to 



= be~P l n(l) 



n{l) 



mbe-P l Y^n{l'). 



(2) 



v>i 



This expression can be solved to generate a recursion 
relation 



n(l+ 1) _ 1 + 1 



el 31 - bl 



I e W+!) - b{l + 1)6 



(3) 



For any population there is a maximum sustainable 
string-length l max , the stability analysis of which 0, 
leads the stationary-state interdependence of b and j3: 



< 



b = 



(4) 
(5) 



For finite-length bit-strings, / max of course cannot exceed 
the bit string itself. 



III. DYNAMICS IN THE SIMPLE PENNA 
MODEL 



Penna model dynamics can be divided into three cases: 
(A) The birth rate or/and mutation rate are altered from 
their steady state values leading to growth or decline of 
the total population n, and associated changes in the dis- 
tribution n{l) of the sub-populations with string-lengths 
I. (B) A change in n(l) from the steady-state distribution 
predicted in equation is followed by relaxation back 
to the steady-state. (C) Within a sub-population n(l), 
the distribution of ages n(x, I) can fluctuate with time if 
the population is not in steady state. 



Steady Growth and Decline in the Simple 
Penna Model 



With an unsuitable choice of birth rate, mutation 
rate and maximum string-length present in a population, 
stationary-state behaviour will not be found. Eventually 
the behaviour of the entire population will be dominated 
by the growth or decline of the longest /-type. The pop- 
ulation can exist in a state of steady growth or decline 
in which the relative sizes of sub-populations remain the 
same. The governing equation can be written as 



(0,0 = be- pi ^nj{x,l) 



(6) 



x=0 



V>1 x=Q 

If we label the rate of growth r, then in any timestep the 
number of young produced is 1 + r times that in the pre- 
vious timestep. Once the population has been growing in 
this manner for some time, populations at successive iter- 
ations are related by %•+! (x, I) = (l+r)rij(x,l). As in the 
steady-state case, n(l) is defined to be the sum over ages 
of n(x,l). We define L r (l) so that nil) = L r il)n(0,l). 
This leads to a simplified steady-growth equation: 







be- 0l rij{l) 



(1 



M0 



■nil) 



(7) 



+ mbe-W^Tnjil') 
v>i 



Where L r (l) is given by 



M0 



l 

l+r 



1 



1 

l+r 



(8) 



The steady growth equation can be manipulated to give a 
recursion relation for the relative sizes of successive rij{l) 



Tljjl + l) Lrjl + 1) 

n,(0 ~ M0 

(l+r)e^ - &M0 
X (l + r ) e W+i) -bL r (l + l)e-P' 



(9) 



The conditions for steady-growth give a limit on the 



value of the maximum sustainable string-length l n 
determine the value of the birth rate b. 
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(10) 
(11) 



In the limit of vanishing r, a power series expansion of 
these expressions will give, to leading order, equations 
(ITO . 
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FIG. 2: Lifespan distributions for a Penna model undergoing 
steady growth with i max = 30, r = 0.05 (+) and i max = 30, 
r = -0.05 (A). 



FIG. 3: Xk (+) and decay time r k (A) plotted against lk for 
a simple Penna model with i max = 30, (3 = 



B. Sub-population Dynamics in the Simple Penna 
Model 

For sub-population dynamics we consider the time-step 
evolution of an arbitrary distribution of n(l) (I < Z max ) 
where the birth and mutation rate take their steady- 
state values. The dynamics of sub-populations within 
the Penna model can be considered to be that of a series 
of decay modes. Any sub-population can be expressed 
as the sum over contributions from a set of decay modes 
such that 

rij(x,l) = }^A k rik,j{x,l), (12) 

k 

where the constants A k are to be determined. Within 
each decay mode, labelled by its index k, the population 
dies away exponentially so that 

nk,j+i (x, I) = (1 - Xk)nk,j(x,l). (13) 

Within each mode there is a maximum value of I above 
which rikj(l) is zero, this is labelled Ik- The mode index 
k is chosen so that k is the number of non-zero sub- 
populations within the mode. Once the relationship be- 
tween rij(x,l) and rij + i(x,l) has been established equa- 
tion Q can be considered to be a sum of eigen-equations 
each of which governs the behaviour of a given mode. 
The equation governing the time-step evolution of the 
longest string-length within a mode can be written as 

oo 

n k , j+ i(0,lk) = be- pik ^n klj (x,l k ). (14) 

The time step evolution of n k j(x, I) is known (|13JI so the 
sum over ages can be evaluated and Xk identified as the 
solution to the equation 



A fc = l- be'? 1 -. (15) 



The characterstic decay time Tk for any mode is defined 
as the time taken for the mode to decay to e" 1 of its 
initial size. 

Manipulation of equation for an individual mode 
gives a recursion relation for the relative sizes of sub- 
populations within a mode 

2^0 + 1) Lk{i + \) 

n kij (l) L k (l) 1 0) 

(I - Xk)eP l - bL k (l) 
X (1 - A fe )e^+i) -bL k (l + l)e-^ 

where L k (l) is given by 



M0 = : _ ■ (17) 

1 l-Xk 

Each mode mimics the behaviour of a population under- 
going steady decline, as given in equation J5Jl. Note that 
the dynamic nature of this model means that the general 
steady-state solution presented in does not give a re- 
cursion relation of exactly this nature. Naturally in the 
limit of vanishing Afe the recursion relation above gives 
that from the steady state Penna model. 

n k j{l) retains its meaning from the steady-state anal- 
ysis as a sum of nk,j(x,l) over ages x. As in the steady 
state analysis, L k (l) gives the sum over ages: n k _j(l) — 
L k (l)n k j(0, 1) The decay mode for which l k = l max has a 
decay rate of (Afc = 0) and is the steady-state solution. 

An arbitrary distribution of rij(l) can be broken down 
into a sum over decay modes. Decomposition of an 
rij(l) distribution into decay modes can be done using 
a 'top down' approach: The largest value of / for which 
rij(l) ^ gives the l k value for the largest mode, A k is 
then chosen so that n,j(l k ) — A k n k j(l k ). Having deter- 
mined A k , A k ^i follows as: rij(l k — 1) = A k n k j(l k — 1) + 
A k _in k _ij(l k — 1). This process is repeated, accounting 
for all contributions from higher modes at each Z, until all 
A k are determined. Any distribution of rij(l) (n(l) = 
for I > Z max ) can be uniquely broken down into decay 
modes. 
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FIG. 4 
1 I 



Sizes of sub-populations within decay modes for /3 = 
:illJ .max = 30 with l k = 25 (+), h = 20 (A) and l k = 15 (x). 
These plots have been rescaled (they are not all normalized 
to 1) to plot them all on the same axes. 



Age Distribution Dynamics in the Simple 
Penna Model 



The dynamic and steady-state behaviour of the sim- 
ple Penna model is dependent on the distribution of ages 
taking a particular form. An arbitrary distribution of 
n(x, I) will have its own dynamic behaviour and must be 
considered as a seperate case. In the simple Penna model 
the evolution of n(x, I) over time without contributions 
from mutation can be dealt with by a Leslie Matrix ^1 
of rank I acting on a vector n(x,l) where the vector com- 
ponents correspond to ages. For n(x, I) the allowed ages 
are 0,1,2,3.../ — 1 and the timestep evolution can be 
descibed by the following matrix equation 



n(l,l) 
n(2,l) 
n(3,0 



V 



n(l,l) 
n(2,l) 
n(3,f) 



/ 



3+1 



V 



(18) 



/ 



where the matrix A is given by 



/ fee-' 3 ' 
1 



v ; 



be- 01 be-? 1 



1 
1 



(19) 



We require eigenvalue solutions to this so that 

n(x,l) j+1 =tn(x,l)j (20) 

£ are the eigenvalues of the matrix A. These eigenvalues 
can be identified as non-unity solutions to the polynomial 



- (be' 131 + l)t l + be- 01 = 0. 



(21) 



The large eigenvalues give decay at roughly the rate pre- 
dicted from the n(l) analysis and oscillation around this 



decay. The small eigenvalues play no significant part over 
large time scales. £ max the largest eigenvalue, is related to 
the decay rate of the decay mode Ukj{l) by Afc = 1— £ roax - 
If / << then the long term behaviour of any pertur- 
bation in n(l) will be dictated by a series of decay modes. 



IV. ANALYSIS OF COMPUTATIONAL 
DYNAMICS 

A Simple Penna model simulation was run with b and (3 
chosen to give Z max = 30. The population was initialised 
with n(x, I) = with the exception of a spike at ro(0, 25) 
where 10 5 organisms were created. The simulation was 
then allowed to run generating data for n(l) at each iter- 
ation. As the maximum lifespan in the population is less 
than Z max the population will eventually decay away to 
nothing. The spike was chosen as the intial distribution 
as any initial distribution can be considered to be a sum 
of such spikes. 

Initially the dynamics of n(25) will be dominated by 
age distribution dynamics and cannot be explained in 
terms of modes. After 100 iterations this noise has all 
but disappeared and the subpopulation can be seen to 
decay exponentially as predicted by equation 113(1 (see 
fig. EJ). All other subpopulations initially grow to later 
enter a period of exponential decay (see figEJ). 

Considering the population after 100 iterations, modes 
with Tfc < 100 will play an insignificant part in the long 
term dynamics of any subpopulation. After large number 
of timesteps t, each sub-population, to a good approxi- 
mation, can be represented by the highest few modes 
(with largest decay times): 

n t {l) w A 25 n25,o(0 e ^ ( 22 ) 

-t 

+ -424«24 : o(0 eT24 

-t 

+ -423"23,o(0 e " 23 

-t 

+ A 2 2n 2 2,o(0 e " 22 ■■■ 

The notation used is n k j(l) where k is the mode index, t 
is the number of timesteps after the 100th iteration. 
are mode coefficients (determined from the distribution 
of n(l) at the 100th iteration). 

Taking more modes into account will give a more accu- 
rate picture so that the behaviour of any subpopulation 
can be described exactly (in the absence of age dynamic 
induced noise) by 



25 



k=l 



CONCLUSION 



(23) 



As a natural extension to our work on the steady-state 
Penna model we have considered and obtained analytical 
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FIG. 5: A plot of n(25) against time from simulation results. 
The population is initialised with n(0, 25) = 10° and all other 
n(x, I) set to zero. The sudden drop in population at 25 
timesteps arises because all of the initial population die once 
they reach age 25. After 100 iterations the noise induced by 
age distribution dynamics has essentially disappeared and the 
population dies away exponentially with time constant r k of 
832.4. 
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FIG. 6: Analytical results for the decay of n(20) and n(10) 
are compared with those from simulation. The upper graphs 
correspond to n(20), the lower to n(10). Analytical results are 
given by the solid lines, simulation results by the dotted lines. 
The analytical results are obtained using only the highest 4 
modes. 



expressions for the various forms of dynamics which can 



be displayed by a simple Penna model. Our approach can 
be applied to the governing equations derived in |£| to 
analyze the dynamics of a variety of Penna models, such 
as those with external death rates and birth cutoffs. For 
multiple disease Penna models (T > 1) the breakdown 
of the population into decay modes remains valid so long 
as there is no bias in the distribution of non-terminal 
mutations within bit-strings. 

The decay coefficients (X k ) and relationship between 
n k ,j(l) within modes for models with T > 1. are given as 
they may be of particular interest. Each X k satisfies 



1 



At- = 1 



l-A, 



1 - 



-be~ p{lk ~ T+1) . 



(24) 



l-A, 



The recursion relation between successive sub- 
populations within a mode is 



n k j(l + l) _ Ci+^Lkil + l) 
n k ,j(l) 



Crp_ 



! L k (l) 

(1 - \ h )eP«- T +» 



(25) 



bL k (l) 



(1 - X k )eP( l + 1 - T + 1 ) - bL k (l + l)e-P' 

where L k (l) is unchanged from the single mutation case 
and given by equation (|17() . 

The dynamic behaviour we have considered is for un- 
regulated populations where there is no Vcrhulst factor 
and steady state is obtained by suitable choice of birth 
rate and mutation rate. Dynamics in the presence of 
a Verhulst factor is considerably more complicated as 
any change in a sub-population will act to change the 
birth (or death) rate, affecting all other sub-populations. 
The modes, which in the absence of a Verhulst factor 
are independent, become coupled through a shared birth 
(or death) rate. The resulting series of coupled, non- 
linear difference equations cannot be generally studied 
using the decay modes analysis we have presented. How- 
ever, depending how the Verhulst factor is implemented 
there may be regimes in which the coupling between de- 
cay modes is sufficiently weak that they may be treated 
as independent. 

The authors would like to thank M. E. Cates and 
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